#include <fstream>
#include <iostream>
#include <cstdlib>
#include <vector>
#include <omp.h>
#include <sys/resource.h>
#include <ctime>
#include <cmath>
#include <algorithm>

using namespace std;

struct Point
{
	double x;
	double y;
	friend ostream & operator << (ostream & ostr, Point & p )
	{
		ostr << p.x <<" "<<p.y ;
		return ostr;
	}
	friend Point operator -  (Point & p1, Point & p2)
	{
		Point p;
		p.x = p1.x - p2.x ;
		p.y = p1.y - p2.y ;
		return p;
	}
	double length()
	{
		return sqrt(x*x + y *y);
	}
	double square_length()
	{
		return (x*x + y * y);
	}
} ;

void printElapsed (char *desc, struct timeval *start, struct timeval *end, int niters )
{
	struct timeval elapsed;
	if (start -> tv_usec > end -> tv_usec)
	{
		end -> tv_usec += 1000000;
		end -> tv_sec--;
	}

	elapsed.tv_usec = end -> tv_usec - start -> tv_usec;
	elapsed.tv_sec = end -> tv_sec - start -> tv_sec;

	printf("\n %s total elapsed time = %1d (usec)",
			 desc, (elapsed.tv_sec*1000000 + elapsed.tv_usec)/niters);

}/* printElapsed() */

double realrand()
{
	return rand () / (double) (RAND_MAX);
}

Point rand_point()
{
	Point p ;
	p.x = realrand();
	p.y = realrand(); 
	return p; 
}

/*The computing domain is bounded by [0, 1]*/
void generate_sites (unsigned * tags, Point * sites)
{
#pragma omp parallel shared (tags, sites)
	{
		int thread_id = omp_get_thread_num();
		for (int i = tags[thread_id] ; i < tags[thread_id + 1]; ++i)
		{
			sites[i] = rand_point();
		}
	}
}

void training(unsigned * tags, Point * sites, 
			  unsigned * weights, unsigned * local_min, 
			  unsigned  num_threads, unsigned num_iterations)
{

	Point seeds ;
	for (int loops = 0; loops < num_iterations; ++loops)
	{
		seeds = rand_point();
#pragma omp parallel shared (tags, sites, weights, seeds)
{
	    // Find minimum on this processor
	    unsigned thread_id = omp_get_thread_num();
		unsigned pos = tags[thread_id] ;
		for (int i = tags[thread_id] + 1; i < tags[thread_id + 1] ; ++i)
		{
			if ((sites[i] - seeds ).square_length() 
				< (sites[pos] - seeds).square_length())
			{
				pos = i; 
			}
		}
		local_min [thread_id] = pos; 
}
        //Find the global minimum 
        //cout<<"SEEDS"<< seeds <<endl;
		unsigned idx = 0;
		for (int i = 1; i < num_threads; ++i)
		{
			//cout<<local_min[i]<<endl;
			if ((sites[local_min[i]] - seeds).square_length() 
				< (sites [local_min[idx]] - seeds).square_length())
			{
				idx = i;
			}
		}
		unsigned global_min_pos = local_min[idx];
		sites[global_min_pos].x = (sites[global_min_pos].x * 
								   weights[global_min_pos] + seeds.x ) 
			/ (weights[global_min_pos] + 1);
		sites[global_min_pos].y = (sites[global_min_pos].y * 
								   weights[global_min_pos] + seeds.y )
			/ (weights[global_min_pos] + 1);
		weights[global_min_pos] ++;
		//cout<<global_min_pos<<endl;
	}
	

}

void divide_group (unsigned size, unsigned num_threads, unsigned * tags)
{

#pragma omp parallel shared(tags)
{

    unsigned group_size = size / num_threads ;

	int thread_id = omp_get_thread_num();
	int d = size % num_threads;
	unsigned start = 0;
	unsigned end = 0;
	if (thread_id < d) 
	{
		start = thread_id * (group_size + 1 ); 
	}
	else
	{
		start = d * (group_size + 1)  + (thread_id - d) * (group_size); 
	}
	tags[thread_id] = start;
	
}
}
int main( int argc, char * argv[])
{

	int num_sites = 0;                  /* num of sites */
	int num_threads = 0;               /* num of threads */
	unsigned num_iterations = 0;

	struct timeval tmstart, tmend;
	struct timeval genstart, genend;
	struct timezone tzp;

	/*----------------------------------------------------------------
	 * Check arguments
	 *---------------------------------------------------------------*/
	if (argc < 3)
	{
		printf("Usage : %s [num of sites] [num of iterations]\n", argv[0]);
		return 0;
	}
	num_threads = omp_get_max_threads();
	if (num_threads > 16)
	{
		printf("Usage : %s [numints]\n - threads overflow. Maxium number of threads is 16.", argv[0]);
		return 0;
	}

	num_sites = atoi(argv[1]);
	num_iterations = atoi (argv[2]);
	printf("\nnum_threads = %d, num_sites = %d, num_iterations = %d\nExcuting %s ...\n", omp_get_max_threads(), num_sites, num_iterations, argv[0]);


	//Step 1. Divide group.
	unsigned  * tags = new unsigned [num_threads+1];
	tags[num_threads] = num_sites;
	divide_group (num_sites, num_threads, tags);
	for (int i = 0;  i < num_threads+ 1; ++i)
	{
		cout<<tags[i]<<' ';
	}
	cout<<endl;


	//Step 2. Choose initial sites ;
	Point * sites = new Point [num_sites];
	generate_sites (tags, sites) ;

	ofstream ostr_ori("ori.in");
	ostr_ori<<"4\n0 0\n1 0\n1 1\n0 1\n"; 
	ostr_ori<<num_sites<<endl;
	for (int i = 0; i < num_sites; ++i)
	{
		ostr_ori<<sites[i]<<endl;
	}
	ostr_ori<<endl;

	//Step 3. Training the sites ;
	unsigned * weights = new unsigned [num_sites];
	for (int i = 0; i < num_sites; ++i ) weights [i] = 1;

	unsigned * local_min = new unsigned [num_threads];
	
	training(tags, sites, weights, local_min, num_threads, num_iterations) ;

	ofstream ostr_ref ("ref.in");
	ostr_ref<<"4\n0 0\n1 0\n1 1\n0 1\n";
	ostr_ref<<num_sites<<endl;
	for (int i = 0; i < num_sites; ++i)
	{
		ostr_ref<<sites[i]<<endl;
	}
	ostr_ref<<endl;
	
	//Step 4. Release all the memories.
	delete []sites ;
	delete []weights ;
	delete []tags ;
	delete []local_min ;
	return 0;
}
